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The collision consistency between the BGK collision model equation and lattice-BGK 
(LBGK) model is proposed by researching the physical significance of the relaxation 
factor r in LBGK model. For microscalar flow in which the continuum hypothesis is not 
still satisfied, the collision consistency r = 1.0 should be ensured when using the LBGK 
model for simulating microflows. The results of simulating microchanncl Poiseuillc flow 
with constant pressure gradient under collision consistency by using LBGK model are 
well consistent with the analytical solutions, and the accuracy of these results is three or 
four orders of magnitude higher than those that don't satisfy the collision consistency. 
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1. Introduction 

Micro-electro-mechanical systems (MEMS) technology has been developed rapidly 
in recent yearsff^ The characteristic length scale of MEMS is typically the order 
of microns, and the ratio of the mean free path to the characteristic dimension 
(i.e. Knudsen number Kn) can not be negligible. The dynamics associated with 
microchannels can thereby exhibit rarehed phenomena and compressibility effects. 
The former is the emergence of a slip velocity at the wall boundary in microchannel 
flows. The latter is the significant nonlinear pressure drop of gas flowing in a long 
microchannel. Because conducting experiments in micrometer-size is a big challenge, 
numerical simulations of MEMS become very important tools of investigation. But 
the methods commonly used in simulation of MEMS, such as molecular dynamics 
(MD), direct simulation Monte-Carlo (DSMcjS and direct numerical simulation of 
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Boltzmann equation, usually requires a tremendous amount of computer time and 
memory. Recently, for its intrinsic kinetic nature, the Lattice Boltzmann method 
(LBE) has been become an attractive method for simulation of m icroscalar flows 
where both microscopic and macroscopic behaviors are coupled 

Form its birth nearly 20 years ago (1988)''', the lattice Boltzmann method 
(LBMjMOl 

has met with significant success for the numerical simulation of a large 
variety of fluid flows, and has emerged as an alternative numerical technique for 
simulating fluid flows. The LBE method usually solves the Bhatnagar-Gross-Krook 
(BGK) model equation, a simplified model Boltzmann equat ion, o n a discrete lat- 
tice, which is usually named the lattice-BGK (LBGK) modelP^IS Si nce th eoretical 
connections between the LBGK model and the Boltzmann equatio n^ ! ^ and the 
preliminary link between LBGK model and the Burnett-type equation!^ were es- 
tablished, the LBE model can be valid for rarefied gas flow provided the Mach 
number is small. Moreover, since the standard BGK equation is able to simulate 
highly non-equilibrium gas flows, the LBGK model should also be applicable to 
rarefied microflows theoretically!^""^ 

One of the most important issues of LBGK model for simulating microchannel 
flows is the introduction of Kn number into LB models. The popular treatment is 
the construction of the relationship between Kn and relaxation time r of LB models. 
Nie et al. used the explicit LBE formulation and related the nondimensional relax- 
ation time r e to Knudsen number Kn as: Kn = a(r e — 0.5)/ pH for a microchannel 
of height H and gas density p. The factor 0.5 attributes to the explicit treatment of 
the collision ter m and a is chose n to best match the simulated mass flow rate with 
experiments P31 Later, Lim et a/ff^ proposed a different relation between Kn and 
r e for the explicit LBE formulation without the correction factor of 0.5 as in Nie 
et al.. Their Kn for a long microchannel was defined as Kn — (6xT e /H) / (po/p), 
where p and po are the local pressure and pressure at the outlet of the microchan- 
nel, respectively. Zhang et al. gave di ffere nt definitions of r e with different constant 
factors among various lattice models!^ Recently, Lee et al. proposed a definition 



of r about fully implicit LBE. 



Another one is the implementation of boundary condition. Nic et al. used the half 
bounce-back rule for the slip effect at the surface. Lim et al. used a specular reflec- 
tion model to generate slip effect. Zhang et al. adopted the Maxwellian scattering 
kernel to address the gas molecule and surface interactions with an accommoda- 
tion coefficient a. Lee et al. proposed a wall equilibrium condition according to the 
assumption of rough surface on the characteristic length of gas molecules. 

For the former issue, the previous research is still related to the framework 
of Navier-Stokes equations. This is because the expression of the kinetic viscosity 
coefficient derived from the LB equation recovering the Navier-Stokes equations, is 
still adopted in the makeup of the relationship of Kn and r. Furthermore, when 
Kn approaches zero, the Boltzmann equation can be reduced to the Navier-Stokes 
equation. But all the existing LBE for microscale flow can not be reduced to the 
Navier-Stokes equation because of using the relationship between Kn and r. 
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The previous researchers always focus on recovering the Navier-Stokes equations 
precisely from the LBGK model and ignoring the deviation of the LBGK model 
and Boltzmann equation. This deviation may not do any effect in simulation of 
macroscale flows. However, the character scale of microflows is much smaller. In 
such a case, keeping the collision consistency of the LBGK model and Boltzmann 
equation should be paid special attention. This collision consistency requires the 
collision frequency of LBGK model is equal to that of BGK model equation. With 
the further study on the relaxation factor r, we found that these two collision 
frequency are equivalent when r = 1.0 and the simulation results are the most 
approximate. This principle is called collision consistency of LBGK model. 

On the other hand, the boundary treatment is based on the Maxwell slip model 
straightforwardly. There is not any adjustable parameter in our model and is similar 
to the Newmann boundary condition in some degree. The implement is not related to 
the Navier-Stokes equations any more, so it can be applied in simulating microflows 
not only in slip regime but also in transition regime. 



2. Lattice-BGK Equation and Collision Consistency 

To begin with, the BGK equation is: 

^ + c-Vh=--(h-h e<l ) (1) 
at r 

here h — h(x,c,t) and h eq = h eq {x, c, t) are the molecular distribution function and 
equilibrium distribution function, respectively. 

D/2 



exp 



m(c — u)^ 
2k B T 



(2) 



where c and u are molecular velocity and macroscopical velocity respectively; tq is 
the collision time; 1/tq represents the collision frequency: 1/tq = c/A. where c is 
molecular mean velocity; A is molecular mean free path. 
LBGK model equation is: 

dfi , _. 1 



dt 
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here fi is particle distribution function: 
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(3) 



(4) 



Ci is lattice discrete velocity and Ui is weighing factor, tq is the collision time. 

LBGK equation ([3]) can be further discretized in space and time. The completely 
discrctizcd form of Eq. ([3]), with the time step At and space step Ax = CiAt, is: 



M x + Cl At,t + At) ~ f t ( x ,t) = -- [/i(x,t) -/f*(x,t)] 



(5) 



where x is a point in the discretized physical space, and it is worth notice that t 
commonly considered as a nondimensional relaxation time, t = To/At, represents 
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the collision consistency of LBGK model by the analysis of the physical signification 
of To and At. 

The above discrete LBGK equation is usually solved in two steps: collision step 
and streaming step. So the physical process in LBGK model can be regarded as 
a large number of particles take place one collision, through the distance Ax and 
the time At. Therefore, Ax and At can be regarded as the mean free path and the 
collision time in the LBGK model, respectively, and 1/At represents the collision 
frequency in the LBGK model. Consequently r = To/At means the ratio of the 
collision frequency of the LBGK model to that of the BGK model. From the gas 
dynamics, the collision frequency to is in direct proportion to the density, and 
in dependence on temperature, but independence of molecular velocity. In order 
to furthest approach to BGK model, the collision frequency of the LBGK model 
should be equal to that of the BGK model, i.e. t = To /At = 1. Especially in the 
microflows (Kn > 10 -3 ), the collision consistency must be much more ensured for 
the collision frequency of molecule in the microflows is far lower than that in the 
macroscopical flows (Kn < 10 -3 ). 

It must be emphasized that LBGK model is between the BGK equation and 
Navier-Stokes equations. We formerly attached importance to recover the LBGK 
to Navier-Stokes equations in the simulation with LBGK model, but now we should 
pay sufficient attention to the consistency of the LBGK model with the BGK equa- 
tion, especially in simulation of the microflows. 

3. Knudsen Number, Lattice Number and Boundary Condition 

In the LBGK model, At is defined as: 

At = (6) 
c 

and in the BGK model, tq is defined as: 



- = \ P) 



Using Eq. © and Eq. ([7]), the relaxation factor t can be determined as: 

_ At _ cAx _ c 1 
tq c A c NhKu 



(8) 



where Njj is the number of the lattice in the characteristic length H, H = AxNh ■ 
Regarding c approximately equal to cP^, Eq. ([8|) can be simplified: 

(9) 



N H Kn 

According the collision consistency, the relation of Kn to the number of the 



lattice is: 
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From Eq. IjlOp. the number of lattice should be in inverse proportion to the Kn, 
rather than be determined randomly. 

Under the collision consistency, Eq. (J5]) can be simplified as: 

/ i (x + c i At,t + At)=/; 9 (x,t) (11) 

The distribution function /j evolved in the lattice actually is the equilibrium distri- 
bution function f^ q . Hence, if the velocity on the boundary can be determined, the 
distribution function of the boundary can also be determined. 

To the slip flow regime, the macroscopical slip- velocity on the boundary can be 
obtained by: 

2- o-, ( du\ 

u s = A — 12 

CT \ d yJ W aii 

where u s is the boundary slip-velocity; a is accommodation factor. Using second 
order difference formula, Eq. (fT2")) . Can be discretized as: 

Kn(4u 1 -u 2 ) , . 

2Ax + dKn 

where u\ and u 2 are the velocity of the first and the second point adjacent the wall, 
respectively. 



4. Numerical Results and Discussion 

In this paper, we study microscalar gas Poiscuillc flows in a constant external pres- 
sure gradient in order to validate the collision consistency. The compressibility effect 
becomes negligible and only the rarefaction effect is accounted for. We study three 
cases with different Kn number value (Kn=0.Q2, 0.05, 0.1), and the number of grids 
adopted in the simulations is chosen according to the collision consistency. Results 
will be compared with the analytical results and the results of more fine grids and 
relaxation factor chosen according to the Eq. |9]). 



4.1. Slip Flows 

In case 1, when Kn=0.l, and the pressure gradient dp/dx — 6.6667 x 10~ 2 . Ac- 
cording to the collision consistency, relaxation factor r = 1.0, and the grid numbers 
are 10 x 150. From the calculating result (see Fig. [I]), we can see the results agree 
closely with the analytical results when t = 1.0, and the relative error is between 
1.43 — 1.54 x 10 -8 . When r — 1/1.4 and r = 0.5, the fine grids are increased, but 
the simulative relative error is very big because of departing away from the collision 
consistency. 

In Figure[2j when Kn — 0.05 and the pressure gradient dp/dx = 6.6667 x 10~ 2 , 
relaxation factor r = 1.0 according to the collision consistency, we can see the 
simulating results agree closely with the analytical results, and the relative error is 
between 1.53 — 1.61 x 10~ 6 . We increased the grids by 1.5 and 2 times and got the 
relaxation factor according to Eq. ([9|), but the relative error was still very large. 



4, 2008 17:4 WSPC/INSTRUCTION FILE tian3 



6 C, Zou et al. 




Fig. 1. The comparison of velocity distribution of micro Poiseuille flow (Kn = 0.1). 



0.012 -, 




Fig. 2. The comparison of velocity distribution of micro Poiseuille flow (Kn = 0.05). 



Figure [3] shows the comparison of calculating results obtained in different grids 
and relaxation factors when Kn = 0.02 and the pressure gradient is 2.0. As can 
be seen from the graph, the results accord closely with the analytical results, but 
in the other two cases, the accuracy of the simulating results decreased instead of 
increasing although the grids increased 1.5 and 2 times, about four times less than 
that of when r = 1.0. 
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Fig. 3. The comparison of velocity distribution of micro Poiseuille flow (Kn = 0.02). 



Table 1. The relaxation factors and relative calculation error in different situations. 



Kn 


Grid 


dp/dx 




Relaxation 
factor 


Relative error 
(Min) 


Relative error 
(Max) 


0.1 


10 x 150 


6.6667 x 10" 


-2 


1.0 


1.43 x 10" 8 


1.54 x 10" 


-8 


0.1 


14 x 210 


6.6667 x 10~ 


-2 


1/1.4 


0.0362 


0.0857 




0.1 


20 x 300 


6.6667 x 10~ 


-2 


0.5 


0.096 


0.102 




0.05 


20 x 300 


6.6667 x 10- 


-2 


1.0 


1.53 x 10" 6 


1.61 x 10- 


-6 


0.05 


30 x 450 


6.6667 x 10- 


-2 


1/1.5 


0.012 


0.05 




0.05 


40 x 600 


6.6667 x 10- 


-2 


0.5 


0.0167 


0.075 




0.02 


50 x 25 


2.0 




1.0 


2.38 x 10" 6 


2.92 x 10- 


(i 


0.02 


75 x 37 


2.027 




1/1.5 


0.0021 


0.02 




0.02 


40 x 600 


2.0 




0.5 


0.003 


0.03 





All the calculating errors of different grids and relaxation factors in these three 
cases are shown in Table 1. For micro Poiseuille flow, the best calculating results 
are obtained since we adopted the relaxation factors in according to the collision 
consistency. So the simulating precision is 3 to 4 times more than that of the other 
results. From Table 1, we can see, the errors are almost the same under the same 
relaxation factors. This also illustrates the influence of the boundary error on the 
whole flow simulating error. 

4.2. Higher Kn Flows 

We now investigate the simulation for higher Kn number in this subsection. The 
ratio of the microchannel length L to height H is chosen as 100, and a 2000x20 
regular grid is applied, which has been verified by grid-dependence. We focus on 
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Fig. 4. Normalized velocity profile at microchannel outlet (_R"n=0.194). 

the streamwise velocity profile U at the exit normalized by the centerline maximum 
velocity and the pressure deviation P' = (p — pi)/p from the linear distribution 
Pi = Po + (Pi — Po) (1 — % I L) . Note that we choose the pressure deviation P' rather 
than the pressure distribution itself in order to show the difference more clearly. 

Fig. [4] shows the velocity profile at the outlet of the microchannel with 
Kn—Q.IQA. Compared with the reliable DSMC result, Nie et al. over-predicted the 
slip velocity obviously and the velocity distribution is more higher in nearly the 
whole outlet cross-section. In our simulation, both the slip velocity and exit ve- 
locity distribution are in consistent with the results of DSMC method. Moreover, 
pressure deviation P' is shown in Fig. [5l The qualitative difference with DSMC 
method has been founded in Nie et al.'s numerical result. Our present result is very 
closely to the analytical solution by Arkilic et al.. And the deviation from the DSMC 
result appears only in the middle part along the microchannels. 

5. Conclusion 

From our study, the physical meaning of the relaxation factor r in LBGK model 
is: the ratio of the collision frequency in lattice system and in the BGK model. 
LBGK model situates between BGK model and macro Navier-Stokes equations 
so the collision consistency needs analyzing when LBGK model is dispersed from 
the BGK collision equation. This is very essential to adopting LBGK model to 
simulate the microflow because the continuum assumption may be break down in 
the microflow. Therefore we should assure the same collision frequency in LBGK 
model and BGK model, that is collision consistency. 
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Fig. 5. Pressure deviation P' along the microchannel (_R"n=0.194). 



Under collision consistency, it is simple to deal with the boundary and determine 
the relation factor. Adopting the collision consistency principle r = 1.0, the results 
of the simulating Poiseuille flow accord closely with the analytical results. Without 
it, the precision of simulating results is three to four times less even if the greater 
grids are adopted. 
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